Preparation

Packages

library(tidyverse) 
library(openxlsx) #read excel file
library(ggrepel) #label text
library(vegan) #Biodiversity Index
library(sf) #Map
library(scales)
library(treemapify)

Call custom function

source("Custom_function.R")

Read data

Use read.xlsx() function to load data from Excel file to R

#library(openxlsx)
filename <- "Plot_NC-BD_20230925_v4.xlsx"
# NuiChua plot
NCplot <- read.xlsx(filename, sheet = "NC_census2")

Display the Structure of data and view some first lines.

dbh: Diametre in mm;
lx,ly: tree coordinates to zero position of quadrat (20x20m);
gx,gy: tree coordinates to zero position of plot (100x100m).

str(NCplot)

## 'data.frame':    17383 obs. of  17 variables:
##  $ PlotID        : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ census        : chr  "C2" "C2" "C2" "C2" ...
##  $ tag           : chr  "01-0048" "01-0049" "01-0049" "01-0050" ...
##  $ stemtag       : chr  "01-0048_0" "01-0049_0" "01-0049_1" "01-0050_0" ...
##  $ spcode        : chr  "dimosp2" "crotth" "crotth" "dryppo" ...
##  $ quadrat       : chr  "0000" "0000" "0000" "0000" ...
##  $ lx            : num  1.119 1.545 1.545 1.188 0.327 ...
##  $ ly            : num  1.54 3.21 3.21 3.44 3.73 ...
##  $ gx            : num  1.119 1.545 1.545 1.188 0.327 ...
##  $ gy            : num  1.54 3.21 3.21 3.44 3.73 ...
##  $ dbh           : num  61 41 28 69 160 ...
##  $ codes         : chr  NA "TNBT" "TNBT" "TBV" ...
##  $ hom           : chr  NA NA NA NA ...
##  $ date          : chr  "2021-03-18" "2021-03-18" "2021-03-18" "2021-03-18" ...
##  $ family        : chr  "Euphorbiaceae" "Euphorbiaceae" "Euphorbiaceae" "Putranjivaceae" ...
##  $ genus         : chr  "Dimorphocalyx" "Croton" "Croton" "Drypetes" ...
##  $ scientificName: chr  "Dimorphocalyx sp.2" "Croton thorelii" "Croton thorelii" "Drypetes poilanei" ...

head(NCplot)

Community (Species composition, size-class structure)

Number of family, genus and species

Use n_distinct() counts the number of distinct value in column: family, genus, scientificName

NCplot |>  
  summarise("family" = n_distinct(family),
            "genus" = n_distinct(genus),
            "species" = n_distinct(scientificName))

##   family genus species
## 1     45    76      93

Number species by family

 (spByFam <- NCplot |> 
     group_by(family) |> 
     summarise(n = n_distinct(spcode)) |> 
     arrange(desc(n)))

## # A tibble: 45 × 2
##    family              n
##    <chr>           <int>
##  1 Euphorbiaceae       7
##  2 Rubiaceae           7
##  3 Phyllanthaceae      6
##  4 Annonaceae          5
##  5 Ebenaceae           5
##  6 Fabaceae            4
##  7 Lamiaceae           4
##  8 Celastraceae        3
##  9 Melastomataceae     3
## 10 Myrtaceae           3
## # ℹ 35 more rows

Figure in treemap layout

library(treemapify)

ggplot(spByFam, aes(area = n,fill = n,label = paste(family, n, sep = "\n"))) +
  geom_treemap(start = "bottomleft") +
  geom_treemap_text(place = "centre") +
  scale_fill_distiller(palette = "Greens", direction = 1) +
  labs(fill = "Species")

Species accumulation curves

Species accumulation curves is a graph recording the cumulative number of species recorded in a particular environment. It used to estimate the number of species in a particular area.

Here, we use b_specaccum() function.

curveNC <-  b_specaccum(NCplot)
head(curveNC)

##   subplot richness    S PlotName
## 1       1 35.24000  400         
## 2       2 48.11000  800         
## 3       3 56.08174 1200         
## 4       4 61.81360 1600         
## 5       5 66.23241 2000         
## 6       6 69.78524 2400

#FIGURE
ggplot(curveNC, aes(x = S, y = richness)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  geom_text(aes(label = round(richness, 1)),
            size = 4,vjust = -0.5, check_overlap = T) +
  theme_bw(16) +
  labs(x = "Survey area (m2)", y = "Number of species")

Forest structure

We usually group trees into 3 size-class of DBH:

small tree [1-5)cm;
medium tree [5-20)cm;
big tree >=20 cm

NCplot <- NCplot |>
  mutate(
    #convert dbh from mm to cm
    dbh_cm = dbh / 10,
    gDBH = case_when(dbh_cm >= 20 ~ ">=20", 
                     dbh_cm >= 5 ~ "[5-20)", 
                     TRUE ~ "[1-5)")
  )

head(NCplot)

We count the number of stems at each diameter class.

(stemByDBH <- NCplot |>
    count(gDBH) |>
    mutate(freq = round(n / sum(n), 3)) |>
    arrange(desc(n)))

##     gDBH     n  freq
## 1  [1-5) 14221 0.818
## 2 [5-20)  2951 0.170
## 3   >=20   211 0.012

Figure

stemByDBH |>
  ggplot(aes(x = gDBH, y = n)) +
  geom_col(aes(fill = n)) +
  geom_text(aes(label = paste0(n, " (", percent(freq), ")")), vjust = -0.5) +
  theme_bw(14) +
  labs(x = "DBH (cm)", y = "Nb of stem", fill = "Nb of stem")

We can also count the number of species at each diameter class.

(
  SpByDBH <- NCplot  |>
    group_by(gDBH) |>
    summarise(
      "family" = n_distinct(family),
      "genus" = n_distinct(genus),
      "species" = n_distinct(scientificName)
    )
)

## # A tibble: 3 × 4
##   gDBH   family genus species
##   <chr>   <int> <int>   <int>
## 1 >=20       26    33      39
## 2 [1-5)      41    71      87
## 3 [5-20)     38    59      72

Figure

SpByDBH |>
  tidyr::pivot_longer(-gDBH, names_to = "rank", values_to = "n") |>
  ggplot(aes(x = gDBH, y = n, fill = rank)) +
  geom_col(position = position_dodge()) +
  geom_text(aes(x = gDBH, y = n, label = n),
            position = position_dodge(width = 0.9),
            vjust = -0.5) +
  theme_bw(15) +
  labs(x = "DBH class", y = "Number", fill = "Rank")

Tree distribution map

Draw a blank plot with draw_blank_plot_1ha() function.

p <- draw_blank_plot_1ha()
p

Add add tree to map

p + geom_point(data = NCplot,
               aes(x = gx,y = gy,size = dbh_cm,color =  gDBH),
               alpha = 0.7) +
  labs(color = "DBH class", size = "DBH(cm)")

For big tree

NCplotA <- NCplot |> filter(gDBH == ">=20")
p + geom_point(data = NCplotA,
              aes(x = gx, y = gy, size = dbh_cm),
              alpha = 0.7,shape = 21,fill = "blue",color = "white") +
  labs(size = "DBH(cm)")

For medium tree

NCplotB <- NCplot |> filter(gDBH == "[5-20)")
p + geom_point(data = NCplotB,
               aes(x = gx, y = gy, size = dbh_cm), 
               alpha = 0.7, shape = 21, fill = "blue", color = "white") +
  labs(size = "DBH(cm)")

For small tree

NCplotC <- NCplot |> filter(gDBH == "[1-5)")
p + geom_point(data = NCplotC,
               aes(x = gx, y = gy, size = dbh_cm),
               alpha = 0.7, shape = 21, fill = "blue", color = "white"
) + labs(size = "DBH(cm)")

For a particular species, px: Dialium cochinchinense

## So do phan bo toan bo cay trong o mau
NCplotSP <- NCplot |> filter(scientificName == "Dialium cochinchinense")

p + geom_point(data = NCplotSP,
               aes(x = gx,y = gy,size = dbh_cm,color =  gDBH),
               alpha = 0.7) +
  labs(color = "DBH class", size = "DBH(cm)")

We can also try with a detail DBH class

NCplot <- NCplot |> 
mutate(
  dbh_cm = dbh / 10, 
  gDBH2 = case_when(
                dbh_cm >= 100 ~ ">=100",
                dbh_cm >= 90 ~ "[90-100)",
                dbh_cm >= 80 ~ "[80-90)",
                dbh_cm >= 70 ~ "[70-80)",
                dbh_cm >= 60 ~ "[60-70)",
                dbh_cm >= 50 ~ "[50-60)",
                dbh_cm >= 40 ~ "[40-50)",
                dbh_cm >= 30 ~ "[30-40)",
                dbh_cm >= 20 ~ "[20-30)",
                dbh_cm >= 10 ~ "[10-20)",
                TRUE ~ "[1-10)")
    )

# Count number of stem in each DBH class
stemByDBH2 <- NCplot |>
  count(gDBH2) |>
  mutate(freq = round(n / sum(n), 3)) |>
  arrange(desc(n))

stemByDBH2

##     gDBH2     n  freq
## 1  [1-10) 16210 0.933
## 2 [10-20)   962 0.055
## 3 [20-30)   162 0.009
## 4 [30-40)    33 0.002
## 5 [40-50)     6 0.000
## 6 [50-60)     6 0.000
## 7 [60-70)     2 0.000
## 8   >=100     1 0.000
## 9 [70-80)     1 0.000

# Figure
stemByDBH2 |>
  ggplot(aes(x = gDBH2, y = n)) +
  geom_col() +
  geom_text(aes(label = n), vjust = -0.5) +
  theme_bw(14) +
  labs(x = "DBHClass (cm)", y = "Number of stem")

Dominant species and Biodiversity Index

Dominant species

The Importance Value Index (IVI) in Ecology is the quantitative measure of how dominant a species is in a given ecosystem. IVI is calculated from 3 values:

relative density of each species (RD)
relative basal area of each species (RBA)
relative frequency of each species (this can be omitted in 1 plot)

\[IVI={\left(relative.density + relative.basal.area \right)}\]

Basal area (BA) of each species can be calculated from diameter measured: $BA = pi * R^2$

Relative basal area (RBA): $RBA = BA / sum(BA)$

Relative density (RD) = number of individuals of each species / total individuals of plot

#   Basal area (BA)
NCplot <- NCplot |> 
  mutate(BA = ((dbh/2)/1000)^2*pi) # unit: m2)

# IVI
IVI <- NCplot |>
  group_by(scientificName) |>
  summarise(n = n(), sBA = sum(BA)) |>
  mutate(
    RBA = sBA / sum(sBA) * 100,
    RD = n / sum(n) * 100,
    IVI = RBA + RD,
    rank = rank(-IVI)) |>
  arrange(desc(IVI))

Figure of IVI

ggplot(IVI, aes(x = rank, y = IVI)) +
  geom_line() +
  geom_point() +
  geom_label_repel(
    data = IVI[1:5, ],
    aes(label = paste0("italic('", scientificName, "')")),
    fill = NA,
    color = "black",
    size = 5,
    label.size = NA,
    min.segment.length = 0,
    segment.linetype = 2,
    segment.size = 0.4,
    parse = TRUE
  ) +
  theme_bw(14) +
  labs(x = "Species rank", y = "IVI (%)")

Diversity Index

A diversity index is a quantitative measure that reflects how many different types (such as species) there are in a dataset (a community). These indices are statistical representations of biodiversity in different aspects (richness, evenness, and dominance).

Richness

Richness is the total number of unique species

Shannon Diversity Index

As species richness and evenness increase, so diversity increases. The Shannon Diversity Index (H) is a way to measure the diversity of species in a community. The higher the value of H, the higher the diversity of species in a particular community. The lower the value of H, the lower the diversity. A value of H = 0 indicates a community that only has one species.

Evenness

The Shannon Equitability Index is a way to measure the evenness of species in a community. The term “evenness” simply refers to how similar the abundances of different species are in the community. This value ranges from 0 to 1 where 1 indicates complete evenness.

\[ Evenness = Shannon Index (H) / ln(richness) \]

Simpson’s Diversity Indices

Simpson’s Index (D) measures the probability that two individuals randomly selected from a sample will belong to the same species

With this index, 0 represents infinite diversity and 1, no diversity. That is, the bigger the value of D, the lower the diversity. This is neither intuitive nor logical, so to get over this problem, D is often subtracted from 1 to give:

Simpson’s Index of Diversity 1 - D

The value of this index also ranges between 0 and 1, but now, the greater the value, the greater the sample diversity. This makes more sense. In this case, the index represents the probability that two individuals randomly selected from a sample will belong to different species.

library(vegan)
# Simpson's Index of Diversity
  D <- diversity(IVI$n, index = "simpson")
# Shannon
  H <- diversity(IVI$n, index = "shannon")
# Richness
  S <- specnumber(IVI$n)
# Evenness
  E <- H/log(S)
# Create data frames 
  res <- data.frame("Richness"=S,"Simpson"=D,"Shannon"=H,"Evenness"=E)
  res

##   Richness   Simpson  Shannon Evenness
## 1       93 0.9014608 2.927733 0.645928

For each DBH class

Create function

The main purpose of creating a user-defined function is to optimize our program, avoid the repetition of the same block of code used for a specific task that is frequently performed in a particular project.

Create calcIVI() to calculator Importance Value Index and figureIVI() to for plotting figure

calcIVI <- function(dataplot) {
  IVI <- dataplot |>
    mutate(G = ((dbh/2) / 1000)^2 * pi) |>
    group_by(scientificName) |>
    summarise(n = n(), sG = sum(G)) |>
    mutate(
      UT = sG / sum(sG) * 100,
      MD = n / sum(n) * 100,
      IVI = UT + MD,
      rank = rank(-IVI)
    ) |>
    arrange(desc(IVI))
  return(IVI)
}

figureIVI <-  function(IVI) {
  ggplot(IVI, aes(x = rank, y = IVI)) +
    geom_line() +
    geom_point() +
    geom_label_repel(
      data = IVI[1:5, ],
      aes(label = paste0("italic('", scientificName, "')")),
      fill = NA,
      color = "black",
      size = 5,
      label.size = NA,
      min.segment.length = 0,
      segment.linetype = 2,
      segment.size = 0.4,
      parse = TRUE
    ) +
    theme_bw(14) +
    labs(x = "Species rank", y = "IVI (%)")
}

Create biodivIndex() to calculator Diversity Index

biodivIndex <- function(IVI) {
  H <- diversity(IVI$n, index = "shannon")
  D <- diversity(IVI$n, index = "simpson")
  S <- specnumber(IVI$n)
  E <- H / log(S)
  res <- data.frame(
    "Richness" = S,
    "Simpson" = D,
    "Shannon" = H,
    "Evenness" = E
  )
  return(res)
}

Run function for each DBH class

Big tree (DBH >= 20 cm)

IVIA <- calcIVI(NCplotA)
head(IVIA)

## # A tibble: 6 × 7
##   scientificName             n    sG    UT    MD   IVI  rank
##   <chr>                  <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Syzygium sp.1             24 3.20  22.8  11.4   34.1     1
## 2 Syzygium sp.2             37 2.05  14.6  17.5   32.1     2
## 3 Drypetes poilanei         24 1.05   7.48 11.4   18.9     3
## 4 Diospyros sp.2            21 1.19   8.44  9.95  18.4     4
## 5 Dialium cochinchinense    18 1.11   7.87  8.53  16.4     5
## 6 Dehaasia sp.              10 0.753  5.36  4.74  10.1     6

figureIVI(IVIA)

biodivIndex(IVIA)

##   Richness   Simpson  Shannon  Evenness
## 1       39 0.9170055 2.919805 0.7969853

Medium tree (DBH in 5-20 cm)

IVIB <- calcIVI(NCplotB)
head(IVIB)

## # A tibble: 6 × 7
##   scientificName        n    sG    UT    MD   IVI  rank
##   <chr>             <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Syzygium sp.2       549 5.64  26.0  18.6  44.6      1
## 2 Drypetes poilanei   521 3.76  17.3  17.7  35.0      2
## 3 Diospyros sp.2      227 1.83   8.43  7.69 16.1      3
## 4 Cleistanthus sp.    327 1.08   4.97 11.1  16.1      4
## 5 Walsura bonii       136 0.868  4.00  4.61  8.61     5
## 6 Diospyros sp.1      116 0.853  3.93  3.93  7.86     6

figureIVI(IVIB)

biodivIndex(IVIB)

##   Richness   Simpson  Shannon  Evenness
## 1       72 0.9062302 2.945215 0.6886708

Small tree (DBH in 1-5 cm)

IVIC <- calcIVI(NCplotC)
head(IVIC)

## # A tibble: 6 × 7
##   scientificName            n    sG    UT    MD   IVI  rank
##   <chr>                 <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Cleistanthus sp.       2767 1.52  22.8  19.5  42.2      1
## 2 Dimorphocalyx sp.2     3215 1.10  16.5  22.6  39.1      2
## 3 Goniothalamus sp.      2230 1.13  17.0  15.7  32.6      3
## 4 Drypetes poilanei       599 0.465  6.97  4.21 11.2      4
## 5 Cleistanthus monoicus   555 0.340  5.09  3.90  8.99     5
## 6 Mallotus sp.            726 0.176  2.63  5.11  7.74     6

figureIVI(IVIC)

biodivIndex(IVIC)

##   Richness   Simpson  Shannon  Evenness
## 1       87 0.8767146 2.706991 0.6061457

Biomass and carbon stock

Biomass

Above-ground biomass (AGB)
Below-ground biomass (BGB)

AGB calculation

Chave et al (2005). Tree allometry and improved estimation of carbon stocks and balance in tropical forests. Oecologia, 145(1): 87–99.

AGB=Total aboveground biomass (kg)
D=DBH (cm)
H=Height (m)
$ρ$=Wood density (g/cm3)
BA=Basal area (cm2)

Habitat	Equation in Chave et al (2005)
Tropical forest DRY with Height	$AGB = exp(-2.187 + 0.916 * ln(ρ * D^2 * H)$
Tropical forest DRY only DBH	$AGB = ρ * exp(-0.667 + 1.784 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)$
Tropical forest MOIST with Height	$AGB = exp(-2.977 + ln(ρ * D^2 * H))$
Tropical forest MOIST only DBH	$AGB = ρ * exp(-1.499 + 2.148 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)$

library(BIOMASS)  

woodensity <- getWoodDensity(
  genus = NCplot$genus,
  species = NCplot$scientificName,
  family = NCplot$family
)

head(woodensity)

NCplot$meanWD <- woodensity$meanWD

Nui Chua plot: Tropical forest DRY only DBH

$AGB = ρ * exp(-0.667 + 1.784 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)$

NCplot <- NCplot %>% 
  mutate(dbh_cm = dbh/10,
         AGB = meanWD * exp(-0.667 + 1.784*log(dbh_cm) + 0.207*log(dbh_cm)^2 - 0.0281*log(dbh_cm)^3))

(AGB <- sum(NCplot$AGB)/1000) #Mg/ha

## [1] 264.7534

BGB calculation

The ratio developed by Mokany et al. (2006). Critical analysis of root : shoot ratios in terrestrial biomes. Global Change Biology, 12: 84-96 offers specific ratios based on forest type and climate zone and are applicable when the aboveground biomass estimate (shoot) is reported at the stand level

(BGB = 0.275 * AGB)

## [1] 72.80718

Total biomass of 1 hecta

(TB <- AGB + BGB)

## [1] 337.5605

Carbon stock

The quantity of carbon can then be estimated by converting biomass to carbon using the IPCC default carbon fraction of 0.47. After that, the carbon stock was multiplied by 44/12 (the carbon atom ratio in the molecular weight of CO2) to get the CO2 equivalent value.

#Carbon stock (ton)
(TC <- TB * 0.47)

## [1] 158.6535

# CO2 (ton)
(CO2 <- TC * 44/22)

## [1] 317.3069

Carbon credit

1 carbon credit = 1 ton of CO2 equivalent (tCO2e) emissions.

In 2023, Vietnam successfully sold 10.3 million forest carbon credits for the first time through the World Bank at a price of $5/ton (Global average selling price ~$11.2)

1 hecta forest at Núi Chúa national park ~ 1586 $

(CO2 * 5)

## [1] 1586.535

Day 01 - Nui Chua Forest Dynamics Plots

Nguyen-Le Xuan Bach

2025-11-30